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We propose an optimized algorithm for the numerical simulation of two-time correlation functions 
by means of stochastic wave functions. As a first application, we investigate the two-time correlation 
function of a nonlinear optical parametric oscillator. 



Q\ ; I. INTRODUCTION 

t-H . 

Spurred by the ever increasing speed of commercially available desk-top computers, the analysis of dissipative 
quantum dynamics in terms of stochastic Schrodinger equation was recently promoted as an interesting alternative 
to the more traditional approach which is based on the solution of the system's master equation. This is because 
workstations are usually strong in in simulating a couple (say 100) runs of a iV-dimensional wave function, but - 
| because of memory limitations — perform rather weak in propagating the 0(N 2 ) matrix elements of the statistical 
operator [Q. In particular in quantum optics, where N may be of the order N = 1000 or larger, stochastic techniques 
enjoy increasing popularity. 

The technique is based on the "unraveling" of a given master equation 
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in terms of a stochastic Schrodinger equation, which propagates a state vector ip r (t) in such a way, that the solution 
of the master equation (0), p(t) — e Lt p(0), is recovered in the stochastic average, 



p(t) = |W(*))<Vv(i)|. (2) 



Here, the subscript r labels a particular realization of the stochastic process, and (. . .) denotes the average over all 
realizations, including a weighted sum over pure initial states. 

A given master equation may be unraveled by a variety of stochastic methods, which involve either appropriate 
Qi generalizations of the continuous Wiener stochastic processes, discontinuous jump processes, or a mixture of both. The 
I* ■ mathematical foundation of the various representations, their relation to a specific physical set-up and their connection 
,£h ! to the quantum measurement problem is nicely reviewed in Refs. [[jjQ. For a semi-popular account, including some 
^ ■ historical and philosophical issues, see Ref. Q. 

Even though the individual trajectories \ip r (t)) enjoy great popularity for the illustration of a quantum systems' 
dynamics, the only physical quantities which may reliably be predicted using the set of \ip r (t)) are the time dependent 

expectation values of quantum mechanical operators, (A(t)) = tr (^Ap(t)j , viz. 



(A(t)) = (Mt)\A\Mt)) ■ (3) 

In particular, multi-time correlation functions may not be computed by means of the representation (|^). The two-time 
correlation function 

(A(t)B(0)) = tr (Ae Lt Bp(0)) , (4) 

for example, can not be evaluated by the recipe "measure B in a suitably selected initial state, propagate the post- 
measurement state to time t, measure A, repeat and average". The Heisenberg operators on the left hand side of 
Eq. (Q) must be treated with special care in open systems with non-unitary time evolution, as algebraic properties 
are generally not preserved by transformation to the Heisenberg picture (i.e. (AB)(t) /= A(t)B(t)) It turns out 
that for non-commuting Ait), B(0), the evaluation of Eq. (||) requires an intricate algorithm involving, at least, one 
pair of coupled stochastic Schrodinger equations. 
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Indeed, the correlator (Q) may be viewed as the single-time expectation value of A with respect to an improper 
state x(t), (A(t)B(p)j = tr LAx(f)), where 

X(t) = e Lt Bp(0) . (5) 

For a pure initial state, p(0) = |^(0))(V'(0)|, the initial value of x(f) is given in terms of a dyadic product of two state 
vectors |-0(O)) and \<p(Q)) = B\if>{0)), 

x(o) = \4>(o))(m\- (6) 

This indicates that the evaluation of (^) may be reduced to the simulation of a single-time expectation value of A with 
respect to x{t), which comes, however, at the expense of dealing with two different wave functions simultaneously. 

In the past, several stochastic wave function algorithms have been proposed which support the numerical simulation 
of a two-time correlation function (HJ^ - . Since all these algorithms are build to yield the correct result in the limit of 
infinitely many runs, they are effective, but they are usually not very efficient. The scheme proposed by Castin et al 
0, for example, is threatened by numerical instability as it relies on the subtraction of possibly large numbers. The 
particular scheme proposed by Gardiner and Zoller H, and Gisin (^], on the other hand, is exponentially inefficient, 
that is the number of trajectories which are needed for a reliable estimate of the desired correlation function is bound 
from below by an increasing exponential function of the correlator's time f. The somewhat more intricate algorithm 
which may be extracted from the work of Breuer and collaborators Q does not suffer from this particular kind of 
ineffeciency, but it is not yet optimized. To date no systematic investigation has addressed the issue of how to tailor 
an algorithm which is both effective and efficient. 

In the present paper we construct an optimal algorithm for a class of stochastic representations which are character- 
ized by pseudo-linear Ito differential equations involving jump processes. The paper is organized as follows. In Sec. || 
we review the unraveling of the master equation for a genuine statistical operator using Ito stochastic calculus. In 



Sec. Ill, we extend our analysis for the stochastic representation of the dynamics of skew-symmetric state operators, 
and we indicate an optimal algorithm. Using the simple example of spontaneous emission of a two-state system we 
analyse the efficiency of alternative algorithm, including the Gardiner-Zoller |j| and Breuer-Kappler-Petruccione 
method. As a non-trivial example, we apply our algorithm to the problem of tunneling in the degenerate optical 
parametric oscillator. The methods of Castin and collaborators @, and of Breuer and collaborators || are revied in 
Appendices |X] and [^. 

II. PROPAGATION OF PROPER STATE OPERATORS 

In this section we seek the unraveling of the model master equation, 

— p — 2apa^ — a^ap — pa^a , (7) 

which - dependent on the algebraic properties of the operators a, - describes spontaneous emiss ion of a 2-level 
atom or the damping of a cavity mode. Generalizations to other systems will be covered in Sec. Ill E| . 
We unravel the master equation (^) in terms of an Ito stochastic differential equation 

\dA(t)) = d£(t)\A(t)) ~ dit(t)tf&\il>r(t)) + dv(t)v\Mt)) , (8) 

where £(f), p(t) and 77(f) are, in general complex, stochastic processes, and differentials are forward oriented, dtp(t) = 
ip(t + dt) — ip(t) etc. Note that the equation is pseudo-linear, as the stochastic increments d£(f), dixit) and drjit) may 
depend on ipr(t)- 

A stochastic process like £(f) is, in general, not differentiable and the stochastic increment d£(f) does not have the 
properties of an ordinary differential, i. e. it can not, in general, be treated as "infinitesimal". However, from Eq. (Q) 
and the assumed smoothness of /5(f) it follows that the stochastic average of the dyadic product |Vv(*))(Vv(t)| must 
be well-behaved and differentiable, 

d (\Mt)) (MQl) = MV>r(f)) (Vv(f)| + WrW) Wr(f)| + Writ)) WM\ ■ (9) 

Note that we are not allowed to drop the last term as would be the case if \dtp r it)) were an ordinary differential. 
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By virtue of Eq. (0) , the differential (||) equals the forward differential of the density operator, dp(t) — p(t+dt) —p(t) . 
For our particular model (Q) 

dp(t) = 2ap(t)a'< dt - a^ap(t) dt - p(t)a^adt . (10) 

Since Eqs. (^) and ( fl0|) must coincide regardless of the value of p(t), they must coincide for any individual member 
of the ensemble p(t), say ip r (t)- Hence we set p(t) = \ip r (t))(ip r {t)\ in Eq. (Q, and in Eq. (^) consider ip r (t) as fixed. 
As we make no assumptions on the action of the operators a and , respectively, we insert the stochastic ansatz (||) 
into Eq. (0) and compare coefficients of the right hand sides of Eqs. (ft) and dlG). The result reads 



(d^t)+d^(t)+dam*(t)) = o 



((de(«) + = dt, ((de(i) + l)d»7*(t)) = 



[dp{t)dp* (t)) = 0, {dfx(t)drj*(t)) = 0, (dr}{t)drf (t)) = 2dt, (11) 

where . . . denotes an average over the stochastic increments for fixed 4>r(t)- 

The unraveling (|^) conserves the norm square |f/v(t)| 2 in stochastic average, as is obvious when taking the trace in 
Eq. (§), 



(A(t)\A(t)) = i ■ (12) 

In fact, all well-established algorithms preserve the norm even for every single realization^]: 

Vt,r: (A(t)\A(t)) = 1 . (13) 

As a rule, algorithms which have been derived from some theory of continuous measurement respect norm conservation 
in the strong version (|l3|), as they deal with proper states which admit a physical interpretation at any time. 

If numerical efficiency rather than merely effectiveness is an issue, algorithms which respect the stronger condi- 
ton (|l^) ought to be preferred. This is because they generate trajectories of equal weigth, avoiding waste of CPU 
time on simulating a large number of trajectories whe n onl y a few "heavy" ones actually contribute to the result. A 
more formal proof of this requirement is given in Sec. 
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The strong conservation law (|13|) imposes additional constraints on the stochastic increments. Taking the stochastic 
differential of Eq. @ we obtain 

(Mt)\dMt)) + (#r(*)|V»(*)> + (<#r(*)l#(*)) = , (14) 

which — using the stochastic ansatz (|J) and dropping the time arguments for readability — amounts to the following 
constraint: 

de + dC - (dp + dp*)Q + drjP + drfP* + d£*d£ - d(*dpQ + dtdrjP 

-dp*d£Q + dp*dpS - dp*dr 1 R + dr]*d^P* - dr]*dpR* + drfdrjQ = , (15) 

where P := (tjj r \a\ip r ) , Q := (ip r \a> a\ip r ) , R :— (ip r \a> aa\il) r ) , and S := (ip r \ a> aa^ a\ip r ) . 

Even if supplemented by the strong condition of norm conservation (|l3|), the conditions ( |TT| ) and ( |l5| ) still offer 
quite some freedom in the choice of a correspondig stochastic model. 

Most prominent are the models of quantum state diffusion jl0| , where the stochastic Schrodinger equation assumes 
the form of a Langevin equation with Wiener stochastic increments, and models of quantum jumps, where the time 
evolution is mostly deterministic, only interrupted by discontinuous jumps at random times. Since jump methods 
seem to be favoured in numerical applications we concentrate on the jump processes in the following. 

In a simple two-branch jump process, at any time t, the tripel (d£, dp, drf) can take on one of two different values, 

(r\P J,, clr,\ - / (£j,MJ>??j) ' with probability dp , , 

{at;,ap, ar)j - ^ ^ ) ; with probab ii ity i _ dp ( w ) 



1 Sometimes this is not obvious at first glance: in some recipies for unraveling ([?]), norm conservation may be violated 
temporarily; however, explicit renormalization will then become necessary, at the latest when expectation values are computed. 
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Most of the time, drj is zero while dt; and d\x take on well-defined, infinitesimal, values, so that the state evolves 
continuously, denoted by the index c. With a certain probability dp oc dt, a "jump" occurs, corresponding to a finite- 
valued tripel fj,j, r]j). Such processes are often derived from a theory of measurement, where jumps correspond to 
"clicks" of some photo detector. 

Note that while Eq. (^) contains stochastic differentials d£, d/i, dr\, we are now left with 6 parameters £j, mj, t?j, 
d£ c , dfi c , and dp, which are deterministic functions of the state vector \ip r (t)}; the only stochastic element in ( |lq ) is 
the decision which branch to take. We have retained the differential notation for the parameters d£ c and dfi Cl which 
will be of order dt, in contrast to £j, /ij and 7/j, which will be of order 1. 

The most simple jump process is obtained by choosing £j = — 1, fj,j = 0; from this, and the strong condition of 
norm conservation (|13|), one readily derives the remaining parameters 

r)j = 1*5- 1 -0^) | — 1 , dp = 2dt |<7|Vv)| 2 , d£ c = y , dfi c = dt . (17) 
The stochastic process (||) now reads 

|# r ) = { (-i + l^hMrV)^), dp (18) 

[ {(^ r \^a\ip r ) - a^a) dt \ijj r ) , 1 - dp , 
and the vector \ip r (t)} will be mapped according to \^ r (t)} i— > |Vv(i + di)), 

r |o- |^r(*)>r X * l^(*)> . dp = 2dt(Mt)\* i *\Mt)) 

mt + dt)) = m {1 rf a u\,i \Mm , i - dp (19) 

This algorithm and its generalizations are quite popular in the quantum optics community, where they have been 
applied to a variety of problems, see Refs. [1-6]. 

III. PROPAGATION OF NON-SYMMETRIC OPERATORS 

A. Construction of an efficient general-purpose algorithm 

As indicated in the introduction, numerical simulation of the two-time correlation function (A(t)B(0)) requires the 
unraveling of the skew-symmetric operator x(t) — e Lt Bp(0). Here, a naive generalization of recipes like ( |l8|) can not 
succeed, as there is no obvious way to define the jump probability for an operator like x(t) which is neither definite 
nor hermitian. 

Recall, however, that the unraveling is only required for a pure initial state, p(0) = \ip(0))(ip(0)\, since the case of 
mixed initial state is easily obtained by means of a suitable weighted sum over such pure states. Denoting 1 0(0)} = 
B\ip(0)}, one has x(0) = \<fi(0))(ip(0)\, which is the dyadic product of two state vectors <f> and ip. Thus unraveling 
may well proceed along the lines of Sec. [ij, this time however for a pair (\4> r (t)) , \ tp r (t))) of vectors, such that the 
skew-symmetric is recovered in the stochastic average: 



x(t) = \Mt)KMt)\, 



(A(t)B(0)) = (Vv(i)|A|<M*)} • (20) 

We start the unraveling of x(t) from an ansatz similar to (j§), but now for a pair of wave functions: 

\d(j> r ) = d£i \(j> r ) — dfiia'a \cj> r ) + dr\\b \4> r ) , 

\dip r ) = d£,2 \ipr) - dfj,2a-'a \ip r ) + d-q 2 a \ip r ) . (21) 
Analogous to Eq. (pd|), we can again derive a set of necessary conditions for the stochastic increments: 



(d£i + d£a+#i<i£a) = 0> ( 22a ) 



(d£i + l)dn* 2 = (d& + l)dm = dt , (22b) 
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d(J,id(j,2 = , 



d/iidr]2 = dn^drji 



drjidrj^ 



0. 

2dt 



(22c) 
(22d) 
(22e) 
(22f) 



We restrict ourselves to a jump process similar to (18), but with free parameters 



( |#r) 
\\d^r) 



I.,: - d/il.cO-tlj) |0 r ) 
J., - d^.cO'tCT) IV'r) 



dp 



1 - dp 



(23) 



As above, the indices c and J denote the continuous branch and the jump branch, respectively. The ansatz ( J23_ 
although less general than in Eq. (pl|), still contains enough degrees of freedom to allow for a wide range of jump 
algorithms, including many of those suggested in the literat ure ; it has the advantage of beeing easy to implement, 

and it will automatically fullfill the conditions (22c) and (22c). 

There are two degrees of freedom in the choice of the free parameters in (^3|) which are mere gauge freedoms with no 
relevance to the efficiency of the process (this has already been observed by Diosi, see ||): for some complex number 
c =/= 0, the transformation 



(d£i,c, d/xi, c , d£ 2 , c , d/j.2,, 



(e*d£i,«c*d/ii iC ,c 1 d&,c,c 1 dn 2 ,c) 



(24) 



will leave (22a) invariant and leads to a process which is equivalent in the sense that it predicts identical trajectories 
for all observable quantities. Similarly, there is another gauge freedom for 771, j, 772, j: for any complex number c ^ 0, 
the transformation 



(»7i,J,f?2,j) 1 — >• (c*1i,J> e 1 J?2,j) 



leads to an equivalent stochastic process. 

We are therefore free to choose the symmetric gauge: 



Vr,t: \Mt)\ = \Mt)\ 



(25) 



(26) 



As this must hold for both the continuous and the jump branch in Eq. (|23|), it fixes the modulus of both gauge 
parameters c and c in J24] ) and Pq). The phase can be fixed by demanding that d£i, 2 ,c and 771, 2, j be real numbers; 
from (22b) it follows that d/ii )2lC , too, will be then real-valued. 

We are now left with 7 real-valued parameters to be determined; since ( |22c| ) and ( |22e| ) are always satisfied by our 
ansatz (p3|), we have to fullfill the following 7 conditions: 



(1 - dp) (d£i,cd£2,c + d£i, c + d&,c) 
(1 - dp)(da,c + l)d/i 2 ,c 
(1 - dp)(d£ 2 ,c + l)d/ii, c 
(1 - dp)(d/ii, c d/i 2 , c ) 
Vijmjdp 
|m,j£|<MI 

I ((d£l,c + 1) - d/^i iC (7 t (j)) \4> r ) I 



(220 



dp (22a) 

dt (H|) 
df 


2dt fl22j|) 

l»»,J*^r>| @ 

|((de 2 ,c + l)-dM2,c^))|V'r)| © 



(22d) 



(27) 



Taking into account that d£i. CJ d^.c, d/Lti iC , d/X2,c and dp are quantities of order dt, while 771^ and 772, j will be quantities 
of order 1, these conditions can be simplified: the values 



dpi, c = dfi 2 



dt 



(28) 



are already determined, and we are left with just 4 constraints for the 5 remaining parameters dp, d£i, CJ d^2, CJ 771J, 

V2J- 
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d£i,G + d&,c = dp (22a) 

mjmjdp = 2dt p 



(29) 

hi,7«r|0r)| = h 2 ,ja|W)| (|2J) 

| ((da,c + 1) - ^Ml.c^)) |0 r )| = | ((rf6,c + 1) - ^M2,c^)) |^r)| @ 

We are thus free to impose exactly one additional arbitrary constraint; since all conditions listed above are either 
required (within our ansatz for the jump process) or mere gauge conditions with no effect on the trajectories of 
observable quantities, it will be this single condition which will determine the efficiency of the algorithm. 

B. Minimizing the error of an jump algorithm 

In a simulation of K trajectories Xj '■= l^j) (V'jli J = 1, . . . , JsT, the true skew operator \ will be approximated by 



\ =iEfc- ( 3 °) 



For any correct algorithm, we have 

X = X- (31) 
At the same time, we would like the error to be minimal: 



X — xl| 2 = minimal . (32) 



If we use the norm ||x|| 2 := tr (x^x), then we find for our ensemble of K trajectories: 



■ ■ 2 

Ix-xll 



K 2 

j, i=i 

For independend trajectories Xj an d Xi, where j ^ I, we have 



tr(x f x), (34) 



where we sum over all realizations r of the stochastic process and denote with P(r) the probability of the realization 
r. Thus, Eq. (|33|) becomes: 



lx-xll 2 




7,1 = 1 / 

j^k / 

-^tr (xrXr) H tr (x f x) - tr (x f x) 

^tr (xtxr) - j^tr (xH) ■ (35) 



The condition (B3) now reads: 



By minimizing 



tr 



(xtxrj = (0r#r) (VvlVO = minimal. (36) 



G 



X;P(r)tr(x&.) (37) 

r 

under the constraints 

^P(r) = 1 and £ P(r) Xr = X , (38) 

r r 

( |36| ) turns out to be equivalent to 

Vr,r': tr(x&0 = tr fe<Xv) • (39) 



We see that it is indeed desirable to avoid algorithms containing rare but "large" trajectories. Unfortunately, ensuring 
Eqs. ( p6| ) or ( |39| ) will be in general impossible, so we replace it by a weaker condition which is necessary but not 
sufficient for Eq. ( |36| ) to hold, but which will nevertheless suffice as the last missing constraint required by the one 
remaining degree of freedom in (x2£ 



V |0) , \$) : J t M4>) (VI V)) = minimal. (40) 



Here, the ensemble average is taken only over the two possible branches in (|23j), not over an ensemble of possible 
states before performing the time step: we want this condition to be fullfilled for every member of the ensemble 
individually. We will show below that this will even guarantee a monotonical decrease of the norm of both \<j>) and 
\ip) in every single realization. 

We can now proceed to derive a jump algorithm fullfilling (p0|). The time derivative can be evaluated separately 
for the jump branch and the continuous branch in (p3|); for the jump branch we find: 

d M4>) (</#)) I a = vlA {<t>\^°\<i>) a>\&*W) ~ (</#> (V#> 

4 ' /i2 $2^2 _ s 2 ^ (41) 



dp 2 
where 

8 := (0|0) = , $ := \*\<t>)\ , * := |*|V>}| ■ (42) 

For the continuous branch, we have: 

d M4>) (v#»i c = 2 (dei lC + d&.c) m) w#> 

-2dt ((4>\^a\cb) (0|0) + {^a\i>) 
= 2 (d£ 1)C + d£ 2 . c ) s 2 - 2dt ($ 2 + * 2 ) s 

= 2dps 2 - 2dt ($ 2 + * 2 ) s . (43) 

So we have on average: 

(eZ«0|0> (VIV))) r = dpd «0|0) + (1 - dp)d ((010) (V#»l c 

= — $ 2 tf 2 + d ps 2 -2di(* 2 + * 2 )« . (44) 
dp v y 

Minimizing this with respect to the only remaining parameter dp yields our choice for the jump probability dp: 

2 

dp = -<5>t>dt . (45) 
s 

The remaining parameters can now readily be derived from (|2 



^i,c = M**^* 2 -±* a W de a , c = - s ^-^ 2 + ^ 2 )dt, ( 46 ) 

d/ii, c = dt , d/^2.c = di 
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C. Comparison with other algorithms 

We are going to compare several alternative algorithms by analyzing the radiative decay of a two-level atom. The 
atom is prepared in its excited state |e) from where it decays into its ground state \g) by means of spontaneous 
emission. The corresponding master equation is given by Eq. (^) with the specification a = \/j\g)(e\. 

We wish to compute the correlation function 

g(t) = (tf(t)&) , (47) 

which for this particular problem should result in g(t) = e -7 *. 
Starting with the initial state 

(|0(O)), |^(0))) =(\g),\e)), (48) 

we observe that the continuous branch in Eq. (^3|) will leave the state unchanged except for scalar factors changing 
the norm of 4> r and ip r . The first jump maps \<fi r {t)) onto 0. As this value can never change again, only trajectories 
with no jump up to time t will contribute to g(t). 



1. The Gardiner-Z oiler algorithm 

In Ref. H an algorithm is suggested which is characterized by the following choice of parameters: 

2 1 
dp = 2* dt, r]i,3 = m,J = -j7 i 

d£i,c = d&,c = ^ 2 dt , dfii, c = dfi 2 ,c = dt , (49) 



with \& = |<t|^>)|, see Eq. (42). Note that this algorithm displays a certain asymmetry in favour of |Vv), as it will 
keep \tp r ) strictly normalized, but exerts little control over the norm of \4> r )- This may pose a serious problem, as the 
following analysis demonstrates. 

Using the parameters (|49"|), trajectories with no jump evolve according to 

(IMt)>,hMt)>) =(|g)e 7t ,|e)) , (50) 

and the probability of getting a trajectory without jumps up to time t is 

P„j = e~ 27t . (51) 

The expectation value 

g{t) = p nJ (VnjWI^I^nj (t)) = (52) 
is predicted correctly, but it requires the simulation of at least 

A^gz » P^, 1 = e 27t (53) 



trajectories to get a statistically significant result. In view of the exponential growth of the right-hand side of Eq. (|53j), 
the algorithm (|49|) is not suited for investigation of the long-time dynamics of g(t). 



2. The doubled- Hilbertspace method 
A different choice of parameters which is symmetric with respect to \ip r ) and \<p r ) is suggested in Ref. Jo) : 

dp = 2dt (('0 r |a t (T|V'r) + (<f> r \^ &\<f>r)) , 7)1,3 = T)2,3 = (( VY | ^ Sr\ 1p r ) + (<j) r \ 6* O | (f> r ) ) , 

C^l.c = d^2,c = ( - 0r|o' t <5'|-0r) + {4> r \<?' a\4>r) , d^l.c = ^M2,c = dt . (54) 

Note that in this formulation the normalization (t/vIVv) + (0r|^r) = lj which is required for t = 0, will be preserved 
by the dynamics. A more detailed account of this algorithm is given in Appendix |A|. 
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For the example of the radiative decay of a 2-level atom, we find that in this algorithm, trajectories with no jumps 
evolve like 



\Mt))\ 1 ( \g) 



(55) 



The probability to simulate such a trajectory is given by 

PnJ^(l + e -^) ■ (56) 

Therefore, one needs 

iVBKP » l + l_ 2lt (57) 

trajectories to predict g(t). 

The analysis of the algorithm proposed in j7|, which we review in Appendix |b| leads to an identical estimate for 
this example. This is no conincidence, as both algorithms are equivalent for any system where <x|</> r ) and cr\ip T ) are 
always orthogonal. 

The estimate (|57|), although already much better than Eq. (p3|), can still be improved, as we shall now demonstrate. 



3. The optimized algorithm 

If we apply our algorithm, defined by the choice of parameters (^5|) and ( ff6| ) , to the example of the spontaneously 
emitting 2-level atom, it turns out that the algorithm degenerates and the evolution becomes deterministic: Trajec- 
tories with no jump evolve like 

\m = e-** | 5 ) , m)) = \e) . (58) 

As dp vanishes for all times, no jumps will ever occur. The trajectory ( |58| ) is the only possible one, and the value of 
the correlation function g(t), 

ff (t) = <*(*)|* + |*(*)> = e~*, (59) 
is predicted exactly by just this single trajectory. 



D. Introduction of a "no-jump-probability" 

As it is desirable not to have to check for jumps in every elementary time interval dt, one needs a more efficient 
way to find the time of the next jump. 

We label the jump times by t\, t%, T3, . . . , and define the quantities qk(t), t > Tk, by 

qk{r k +t) := P( "no jump in }r k ,T k + t]" ) , (60) 

where P denotes a probability. Clearly, the qk(t) are monotonically decreasing, starting from qk{Tk) — 1, and they 
obey the differential equation 

q k {t + dt) = q k {t) {1 - dp(t)) . (61) 

For the time Tk+i of the next jump after the jump at Tk, we have for t > Tk, 

dP(T k+1 £]t,t + dt]) =q k {t)dp(t) = \dq k (t)\, (62) 

(Probability of the jump not occuring before t, multiplied by the probability of any jump occuring in dt). If r is a 
random number with uniform distribution in [0, 1[, then we find that for any t > Tk, 

dP{re [qk{t),q k {t)+dq k {t)[) = \dq k (t)\ = dP (r fe+1 e]t, t + dt]) . (63) 

The inverse function (r) can therefore be used to map uniformely distributed random numbers r G [0, 1[ onto 
jump times Tk+i with the correct distribution. This allows to simulate the process ( p3| ) by integrating equation ( |6~l| ) 
in parallel to the continuous branch of (p3|), and performing a jump whenever q k {t) drops below a previously chosen 
random number r € [0, 1[. 
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E. Generalization for arbitrary master equations 



*(*) = E ( 2& kx(t)at - *ta k m - m4*k) - * [h, m 



(64) 



If we generalize this to a master equation containing several loss channels, as well as an Hamiltonian interaction, 
d 

x(t) = 

fc=i 

with 

X(0) - |^(0)) (-0(0)1 , (65) 

we arrive at the following recipe for unraveling: 
Initialization: 

• Use a gauge transformation |0(O)) h-> c|0(O)), |V>(0)) h-> c" 1 |^(0)}), such that (0(0) 1 0(0)) = (<0(O)|V>(O)); 

• initialize the real-valued auxiliary variable q(0) := 1; 

• choose a random number r G [0, 1[. 

The continuous evolution between jumps is governed by the following pseudo-linear equations of motion: 

2 
s 





d 






d 


\<f>r) 


dt 


d 


|VV> 


~dl 



la k + iH |0,,) 



where, as above, 



— ?V$A, 

k=l 

k=l V ~ / \fc = i 

- s £ - \& k + \A) - [E + i^J K 

:= (0 r |0r) = (lpr\lpr) > $k ■= \&k |0r)| , *fe := |<5fc|</v)| . 



(66) 



(67) 



Whenever the variable q drops below the previously chosen random number r, a jump occurs according to the following 
scheme: 

• Choose randomly one jump operator a k using the statistical weights $fc^fc for the individual channels; 

• apply the following map: 



q< — >1, |0r) 
• choose a new random number re [0, 1[ . 



-a k |0 r ) , |Vv) 



-O-fc \4> r ) 



(68) 



F. Some properties of the proposed algorithm 

From the map ( |6^ ) it can be seen that jumps have no influence on the norms of \4> r {t)) and \ip r (t)) (which is the 
local equivalent to the global condition ([39])). The continuous evolution, however, will cause the norms to decrease 
montonically: 



d n ( 1 11 \ 

- (0 r |0 r ) = 2 (0 r | J2 + ~ - **** J ^ 



k=l 
( 



^E 

fe=l 



y ^0 r |<7fc<7fc|0 r ^) ^Vv|<5fc0'fc|VV 

-| ^0 r | <7£<7 fc |0 r \ + (ll> r \ala k \i>r}) 



(69) 
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From the properties of geometric and arithmetic means, we find: 

j^r^r) = j t (A\A) < 0. (70) 

According to Eq. (|35|), the absolute error is bounded from above by K~ x (cf) r \(f> r ) (ip r \ip r ) , and from Eq. (|70| ) it follows 
that this upper bound is decreasing with time. 

An implementation of our algorithm is available via web-download, see Ref. p| . In order to demonstrate its 
reliability, we have simulated the two-time correlation function of a driven 2-levcl atom. The coupling to the classical 
driving field is described by the Hamiltonian H = ^QQe)(g\ + \g)(e\). Spontaneous emission is described by the 
Lindblad operator (^), with a — ^/7|g)(e|. In Fig. |we compare the results of our simulation with the exact solution. 
We find excellent agreement between the results of the simulation and the exact formula for both the correlation 
function. A somewhat larger model system is considered in the next section. 



IV. THE OPTICAL PARAMETRIC OSCILLATOR 



We consider a system of two resonant optical modes h\ (fundamental) and 0,2 (second harmonic), interacting in a 
X^ 2 )-medium via the Hamiltonian 



0,2 



where k is the coupling constant. The second harmonic mode is pumped by a coherent field of amplitude e: 

Hp = i ^eaJ, — f*ai 

and both modes are damped, with rates 71 and 72, respectively: 



Cp = 71 ^2aipa\ — a\aip — pa\aij + 72 {2a2pa\ — a 



t& 2 p- 



pa} 2 a,2 



(71) 



(72) 



(73) 



A detailed discussion of this system can be found in Ref. @. For pump field amplitudes above the threshold 
e t h = 71 72/^, there exist two classical steady state solutions for the fundamental mode amplitude: 



ai = ± 



;(e- eth) 



(74) 



The quantum state resembles this classical solution: above threshold, the Wigner function consists of two peaks close 
to the classical solutions. The superposition of these peaks is incoherent, so the steady state is a classical mixture 
of two states both of which are well localized and have a well-defined quantum phase. However, tunneling between 
the two stable states is possible, corresponding to a change of the phase angle by n. Therefore, above threshold the 

tunneling events will be the main reason for the decay of the correlation function g(t) := ^ a , 1 }*} a ^ in the long-time 

regime. Following |l3[ ], the tunneling is modeled to be a telegraph noise process with a characteristic time T, which 
leads to an exponential decay 



(75) 



By adiabatic elimination of the second harmonic mode (the pump mode), which is a g ood approximation only in the 
case of a rapidly decaying pump mode, i. e. for 72 ^71, Kinsler and Drummond |13|] found the following analytical 
expression for the tunneling time T: 



T 7F 



A + (7 



71 V A(A - a) 



aj 2 , G 2 



A a ■ a In I — 

(j 



(76) 



where A := e/eth is the normalized pump field amplitude, G := ^^^7172 is the scaled coupling constant, and 
(j := 1 — G 2 /2. Eq. ( [76] ) is obtained in the "potential-barrier approximation" which is valid in the limit of large 
threshold photon numbers, i. e., for G << 1, and a large potential barrier, i. e., it fails for A very close to or below 
the threshold A = 1. 
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We carried out numerical simulations of the full quantum dynamics specified in Eqs. (|7l|-|73|), i- e. without adiabatic 
elimination of the pump mode, using a truncated number state representation of the state vectors. Classical analysis 
predicts a photon number N2 = 1 for the pump mode above threshold, and N\ — 8 for the fundamental mode at A = 2 
for the choosen values 71 — k, 72 = 4k. Since photon number fluctuations are enhanced by the nonlinear interaction, 
the truncation values of a number state representation have to be chosen considerably higher than these values. We 
used truncation values for the fundamental mode between a photon number Ni — 24 for Awl, and Ni = 48 for 
A ~ 2. Tests with truncation values up to Ni = 64 indicated that no significant error was caused by this truncation. 
The Hilbert space of the second harmonic mode was truncated at a photon number N% = 16. 

Results of our numerical simulations of g(t) are depicted in Fig. ^ and Fig. ||. Quite generally, g(t) displays a fast 
transient, which decays on a small time scale oc Key^ , followed by a slow exponential decay oc exp(— 2t/T) which 
governs the bevaviour of g(t) in the limit t — ► 00. 

In Fig. || we depict the tunneling times T which we extracted from out numerical data as a function of the 
normalized pump field amplitude A. We see a good agreement between our simulation and the prediction of Eq. (|7q), 
for intermediate values of A. For larger values of A, the tunneling times predicted by our simulations are consistently 
shorter than predicted by Eq. (|7^) . This may be due to the full quantum dynamics of the pump mode which is taken 
into account in our simulations. 



V. SUMMARY 



We have derived an efficient method for the numerical simulation of quantum mechanical two-time correlation 
functions which is based on stochastic wave function propagation. In comparison with other algorithms, our algorithm 
was demonstrated to be generally more efficient, i. e. requiring less runs for a reliable prediction of g(t). We have 
successfully applied our method for the simulation of g(t) of a nonlinear optical parametric oscillator. 

The tests indicate that our algorithm is stable and efficient, even for "large" (i. e.jhaving large Hilbert spaces) 
problems: despite the seemingly more complicate pseudolinear equations of motion (p6|), a proper implementation 
need not be slower but can be of of faster convergence than previously published algorithms. This stems from the 
fact that the most time consuming operations are the application of operators to vectors (operations of the type 
"matrix times vector"), and all algorithms require the computation of cr^crfe |$) and <j\&k for every computation 
of a time derivative to be used in the numerical integration of the continuous (no-jump) part of the evolution. All 
other quantities required in our algorithm can be obtained by computing scalar products of vectors and purely scalar 
operations, which are comparatively cheap operations. 
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APPENDIX A: THE DOUBLED-HILBERTSPACE METHOD 

The algorithm proposed in M is based on the observation, that the problem of propagating the skew-symmetric 



i((t), as defined in Sec. [II A, can be reduced to solving a master equation for a positive symmetric operator in a 
Hilbert space of twice the dimension of the original space. 

The new density operator £l(t) in this doubled space is defined by the initial condition 

firm- fW°)>M°)l l<M0)>W0)|\ , An 



d Cl(t) = CCl(t) :=( L a J] m- (A2) 



and the equation of motion 

The superoperator C is again a Lindblad operator: 

C Cl = 2£f2£ t - E+Efi - OE+E , (A3) 
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where 

i i) ■ ^ 

The initial value Q(0) can be written as a symmetric dyadic product, 

n(0) = |T(0))(T(0)| , (A5) 

where 



\T(0))-( im) ) 



(A6) 



is a normalized pure state in the doubled Hilbert space. Standard algorithms can now be used to unravel ( |A2| ) into a 
stochastic process for a state vector \T r (t)} in this doubled Hilbert space, e. g. the jump algorithm ([T§| ) can be used: 



| (/Ti( , V) = ( , ^-l+E|T r (*)) Sj|T r (t)), dp = 2dt(? r (t)\^X\r r (ty / (A7) 
((T r (t)|Sts|T r (t)) - Ets) A |T r (t)) , 1 - dp. 

The jump probability is computed as the arithmetic mean, 



dp = 2<it(J r (t)|£ t £|T r (i) / 

= 2dt ((^(t)|ata|^ r (t)) + (Mt)\&&\Mt))) > (A8) 



as opposed to the geometric mean suggested in our algorithm (|15|). However, unlike (|49|), this algorithm is symmetric 
with respect to the two components \<p r (t)) and \tp r (t)). Moreover, the condition 

(Mt)\M*)) + (Mt)\M*)) = i> (A9) 

which is the proper normalization of vectors in the doubled Hilbert space, will be fullfilled for all realizations, so 
exponential growth of the nor m o f one of the vectors as in (|50|) cannot occurn. 



Noteworthy, the condition (A9) is actually more restrictive than necessary: in the doubled Hilbert space, we are 



only interested in expectation values of operators of the Form 



A = ( ° °] : .Ann 



in particular, the expectation value of the identity operator 



T = (J D ) ,A111 



is irrelevant, and norm conservation in the stochastic average, 



(Mt)\Mt)) + (Mt)\Mt)) = 1 . ( A12 ) 

is not required. Therefore, condition ( |A9[ ) is less well justified than the corresponding condition (|l^) in the unraveling 
of master equation (fy for a symmetric operator, and we believe that it should be replaced by condition (f4(i|). 



2 Consequently, a constant factor must usually be applied when computing expectation values, to compensate for this 
normalization. 
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APPENDIX B: THE M0LMER-CASTIN-DALIBARD ALGORITHM 



Already in Q, the following stochastic jump process was proposed (although formulated as in Eq. (B2)) 

I'M*)) 



\\dMt))J 



-l+|CT|^r(*)+^r(t)>| 



\Mt)} 

\Mt) 
\Mt)) 



dp = 2dt \&\ip r (t) +f<M*))T 
, 1 - dp 



(Bl) 



Here, v is a phase factor with \v\ = 1, £ is the operator defined in Eq. (A4), and we use the abbreviation 
\4> r (t) + vip r {t)) = \<j) r (t)) + v\ip r (t)). The vectors are normalized such that \\<j> r (t) + vil> r (t))\ = 1; it is easy to 
verify that the process (Bl) will conserve this normalization. In the original formulation of this method fij, it was 
suggested not to propagate the pair (\ip r (t)) , \<fi r (t)}), but the linear combination \ip v ,r(t)) '■= \<t>r{t) + vip r (t)); this is 
possible, as all coefficients appearing in both branches of Eq. (Bl), as well as the jump probability dp, are functionals 
of \ipi/,r(t))i and we can write: 



-1 



|, /( ,, r(/ ,\ > y - 1 l°"IVv(*))l 1 o-)|^,r(i)) I dp = 2dt\a\4>^ r {t))\' 

\ ((^ v ,r(t)\^a\i} Ujr (t))-aU)dt\^ I/jr (t)) ) l-dp 



(B2) 



However, merging of both vectors \ip r ) and \<f> r ) into one only allows for a simulation of 



Xu{t) := |^, r (t)) (^,r(*)| 



|^ r (t)) (Vr(*)| + V\M*)) (M*)\ + V*\Mt)) (Mt)\ + I0r(*)> (0r(t)| 



(B3) 



To extract the quantity of interest, x(t) = \<p r (t)) (ip r (t)\, one has to perform the simulation for four different values 
of v. 



1 



(x+i(t) - X-i(t) - iX+i(t) + iX-i(t)) ■ 



(B4) 



Since (B2) is just the simple jump process (Jig), applied to \ipt/,r(t)), the problem of propagating the skew object x(t) 
has been reduced to the problem of unraveling 4 proper, positive definite density operators. 

Clearly, the extra work can be avoided (at the cost of having to pro pagate two separate vectors instead of a single 
one) by not merging the vectors and using the stochastic process (Bl) directly. 



APPENDIX C: A CONCIEVABLE IMPROVEMENT WHICH TURNS OUT TO FAIL 



The algorithm derived in Sec. Ill A was based on the requirement to simulate all operators of the type ( |A1C| ), 
disregarding more general operators in the doubled Hilbert space. 
If all we want to compute is the expectation value 



(t) = U r (tM\Mt) 



(Cl) 



for just one special operator A, we could ask whether one can formulate an even more specialized algorithm fitted to 
exactly this problem. More specifically, instead of the error ( [32] ) of the estimate of the skew operator \, we could wish 
that the mean square of the error of our estimate of g(t) should be minimized. Again, we cannot obtain the global 
minimum and have to resort to a local condition, analogous to Eq. (EQ): 



minimal . 



(C2) 



An algorithm fullfilling (C2) can be derived in much the same way as from Eq. (|40|). In particular, the corresponding 
jump probability is: 
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dp = 2]T 



k=l 



ip r \alAa k \(f) r 



1p r \A\(f> r 



dt. 



(C3) 



Jumps will not change the value of |<7 r (f)| in this algorithm^. At first glance, this method looks promising: when 
applied to the simulation of the correlation function g(t) — (o^ct) of a 2-level atom, then we have to choose A = and 
(C3) will contain a product oo = 0. Therefore, no jumps will ever occur for any states \<f> r ), |Vv)i and the algorithm 
becomes deterministic even for the driven 2-level atom. Unfortunately, the results will be incorrect. In general, the 



method derived from condition (C2) turns out to be highly unstable at best, and produces completely wrong results in 
some cases. This happens because it does not posses a property analogous to J70|); rather, |<7 r (£)| rnay very well grow 
large for certain trajectories during the deterministic continuous evolution. Moreover, the algorithm will completely 
avoid to jump into states for which g r (t) vanishes, although such states can be important as g r {t) may very well take 
on finite values at a later time. 



3 Even the phase of g r (t) can be conserved if we relax the requirement of our parameters 771, j and 772,.: to be real-valued. 
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FIG. 1. Time evolution of the correlation function g(t) = {cr* (t)a) for the driven 2-level atom, with Q, — 87. The inset 
shows the corresponding spectrum. The solid line is the result of a Monte Carlo simulation with 5000 trajectories, while the 
crosses represent exact analytical values. 
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FIG. 2. The correlation function g(t) of the optical parametric oscillator for system parameters 71 = k = 1, 72 = 4k, and 
a normalized pump amplitude A = = 2.0, derived from 500 simulated trajectories. The inset displays a close-up view for 
different values of A, showing the fast initial decay before the slow tunneling regime. 
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FIG. 3. Comparison of numerical values for the tunneling times of the optical parametric oscillator and the analytical 
approximation [T6l for system parameters as in Fig. 0. 
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